Dynamic exchange-correlation potentials for the electron gas in 
dimensionality D = 3 and D = 2 
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Recent progress in the formulation of a fully dynamical local approximation to time-dependent 
Density Functional Theory appeals to the longitudinal and transverse components of the exchange 
and correlation kernel in the linear current-density response of the homogeneous fluid at long wave- 
length. Both components are evaluated for the electron gas in dimensionality D — 3 and D — 2 
by an approximate decoupling in the equation of motion for the current density, which accounts 
for processes of excitation of two electron-hole pairs. Each pair is treated in the random phase 
approximation, but the role of exchange and correlation is also examined; in addition, final-state ex- 
change processes are included phenomenologically so as to satisfy the exactly known high-frequency 
behaviours of the kernel. The transverse and longitudinal spectra involve the same decay channels 
and are similar in shape. A two-plasmon threshold in the spectrum for two-pair excitations in D = 3 
leads to a sharp minimum in the real part of the exchange and correlation kernel at twice the plasma 
frequency. In D = 2 the same mechanism leads to a broad spectral peak and to a broad minimum 
in the real part of the kernel, as a consequence of the dispersion law of the plasmon vanishing at 
long wavelength. The numerical results have been fitted to simple analytic functions. 
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I. INTRODUCTION 

The plasmon dispersion relation and the dynamic 
structure factor of the conduction electrons in simple 
metals are known from Electron Energy Loss and In- 
elastic X-ray Scattering experiments.ErB Both electron- 
gas correlations and band structure enter in determin- 
ing these properties and Time-Dependent (TD) Density 
Functional Theory (DFT) provides a general framework 
which can account for both. The generality of the DFT 
methodHu in dealing with inhomogeneous electron sys- 
tems motivates interest in the derivation of sensible ap- 
proximations to the dynamic exchange and.-correlation 
(xc) potential, beyond the adiabatic regimdJia whose ap- 
plicability is limited to low-frequency phenomena. 

A local density approximation for a scalar xc poten- 
tial in time-dependent phenomena was proposed in early 
work of Gross and KohnBllj However, detailed analysis of 
the constraints coming from basic conservation laws has 
shown that inconsistencies can arise and are associated 
with the non-existence of a gradient expansion for the 
frequetLG¥-jdependent xc potential in terrns-|G>£,the density 
alone. □'tlrtil Recently Vignale and KohrJi3t3 have over- 
come these difficulties by resorting to a dynamic xc vector 
potential even in the case of an inhomogeneous system 
subject to an external scalar potential. They obtained an 
explicit local-density expression for the xc vector poten- 
tial in the linear response regime in terms of correlations 
of longitudinal and transverse currents in the homoge- 
neous electron gas. This expression becomes exact when 
the equilibrium electron density and the external poten- 
tial are slowly varying in space, on length scales set by 
kp 1 and vf /oj where kp and vf are the local Fermi wave- 
number and velocity. 



The results of Vignale and Kohn have been inter- 
preted in terms of complex, frequency-dependent vis- 
coelastic coefficients, allowing a non-linear generalization 
up to second order in the spatial .gjjadients.EZI Within 
this framework, Ullrich and Vignaldla have developed a 
simple computational scheme for the calculation of the 
linewidths of non-Landau-damped collective excitations, 
whose decay is solely due to dynamical xc effects. The 
order-of-magnitude enhancement over the homogeneous- 
gas approximation reported in the application to quan- 
tum stripalS emphasizes the relevance of the interplay 
between inhomogeneity and dynamical xc effects. 

The work of Vignale and Kohn has brought new inter- 
est to the evaluation of the longitudinal and transverse 
components of the dynamic xc kernel [f^ c (uj) and f^ c (uj), 
say] in the homogeneous electron gas at long wave- 
length. An interpolation between the known asymptotic 
behaviours at low and high frequency was proposed for 
the longitudinal component by Gross and KohnBE3 and. 
later extended to finite wave, number by Dabrowski.LS 
It was subsequently noticedEHl that the long-wavelength 
longitudinal spectrum is closely related to processes of 
excitation of two correlated electron-hole pairs in the 
electron gas. These have been studied— by perturbative 
methods in the-eady work of DuboisE!rE3 and by several 
other authors .circa The inclusion of dynamic screening 
in the Random Phase Approximation-/RPA) leads to a 
mode-coupling form of the spectrumRa and inclusion of 
final-state exchange processes is needed to recover the 
exact rWgh- frequency behaviour calculated by Glick and 
Long.ca A similar expression has been more recently ob- 
tained by Neilson et a/Ej within a memory function for- 
malism for the electron gas in dimensionality D = 2. 

With regard to the transverse component of the long- 
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wavelength xc kernel fj c (uj), Vignale and Kohnli§£3 have 
given a first-order perturbative estimate and obtained 
the high-frequency limit. We have subsequently briefly 
reported on the results of a full calculation of d^ c {uj) for 
the electron gas in three spatial dimensions,^ _and on 
preliminary results in the two-dimensional case.Lj These 
calculations were based on the two-pair model treated in 
the RPA and corrected for final-state exchange processes 
(hereafter indicated as RPAE). 

In the present work we give a full account of our ap- 
proach to the evaluation of the dynamic xc kernels f^ c (uj) 
and /J c (w), including an exact expression for /J^! (&>) in 
terms of four-point response functions. We also extend 
our calculations to (i) fully evaluate the dynamic xc ker- 
nels for the electron gas with e 2 /r interactions in D = 2, 
and (ii) examine the role of including exchange and cor- 
relation in the screening processes entering the RPAE. 

The lay-out of the paper is as follows. Section |l| 
presents the theory underlying our calculations, with the 
help of three Appendices. We start from the definition 
of fxc T (uj) in terms of the current-current response func- 
tions for the ideal and the real electron gas and proceed 
to evaluate the ideal-gas response at high frequency and 
to derive an exact expression for the real-gas response 
in terms of a four-point response function. An approxi- 
mate decoupling of this latter function into products of 
two-point response functions introduces the two-pair ap- 
proximation. Screening at the RPA level or better and 
phenomenological inclusion of final-state exchange finally 
lead to the formulae used in our calculations. The numer- 



ical results are presented in Section III together with fits 
to analytic functions incorporating the known asymptotic 
behaviours and aimed at facilitating numerical applica- 
tions within the time-dependent DFT formalism. The 
role of exchange and corre lati on in the treatment of each 
pair is studied in Section IV. We conclude with a brief 
summary in Section 



II. THEORY: EXACT RESULTS AND TWO-PAIR 
MODEL 

The longitudinal (L) and transverse (T) kernels 
fxc fa) °f the homogeneous electron gas are defined as 
the k — > limit of the functions 



/ x ^(fc, w ) = — 



X° L , T (k,uj) +n/m 



XL,T(k,ui) + n/m 



L.T 



(1) 



Here, Xl(t) is the longitudinal (transverse) current- 
current response function of the homogeneous fluid at 
density n, X°l(t) 1S tne corresponding ideal-gas response 
function, v% is the Coulomb potential (v% = A-Ke 2 /k 2 



in D 



3 and v% 



2ire 2 /k in D = 2) and v 



I = 0. We 



remarkllS that xl {k, ui) is related to the density-density 
response function x(k,u>) by 



Xh{k,uj) + 



k 2 



(2) 



and that f£ c is proportional to the local field factor 
G(k,co) entering the dielectric function of the electron 
gas, according to f£ c (k,u)) — —v^G{k,oS). 

The longitudinal response function x°(fc,w) is-iranie- 
diately obtained from the Lindhard susceptibilityEjLj in 
D = 3 and from the Stern susceptibility!^ in D = 2 
by using (||) for the ideal Fermi gas. The calculation of 
XTp{k,ui) for the ideal Fermi gas in D = 3 and D = 2 is 
reported in Appendix 

Equation ([!]) may be written in terms of the proper 
current response functions XL.T{k^) 1 



fh T (k,u>) 



1 



X° )T (fc,w) +n/m 
1 



XL,T(k,u) + n/m 



(3) 



This emphasizes that the plasmon does not contribute 
to f^ c (k,u>), leaving only contributions from multi-pair 
excitations in the limit k — > 0. In fact, to leading order 
in the long wavelength limit we have from (p|) 



Im/i^H = Jim 



TO 2 W 2 



k^o n 2 k 2 



Imxi,r(fc,w) 



(4) 



We proceed below to evaluate the imaginary part of the 
kernels, from which their real part will be obtained by 
means of the Kramers-Kronig relation. The equation of 
motion for XijO*-, w) can be obtained from the definition 
of a general response function in terms of unequal-time 
commutators, 



* / e l{u+le)t {[A{t),B{Q)])dt (5) 



U;B)) u 



where A(t) — e lHt Ae~ lHt , (...) denotes a ground-state 
expectation value, and e is a positive infinitesimal. The 
current-current response Xijfe, = {{j^jsL^ui satis- 
fies the equation of motion 



w 2 Xy(k,w) = ( [jl,H],j 



.\ k .ii.r k .ii 



(6) 



where the first term is real and independent of w. It 
will be evaluated below in the discussion of the real part 
of f xc (Eqs. ^5[]l7| ). Quite lengthy calculations, which 
are briefly reported in Appendix [BL lead to the long- 
wavelength result 



Imxij (k, u) 



1 



m 2 V 2 u>' 



Y. Im ((jqP-q;jq'/ 5 -q' 



x r a (q,k)r^(q',-k)+ (fc 2 ) 



(7) 
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where the coefficients r l '(q, k) are defined in Eq. (B7) 
and summation over repeated indices is understood. This 
expression, together with (0), gives an exact result for 
Im./£ T M. 

We now discuss some approximate evaluations of 
Eq. (0). Use of the ideal-gas four-point response func- 
tion in the RHS gives the exact second-order perturba- 
tive value for /.^(w), which is expected to be accurate 
at high frequency (see App. ^). In Appendix ^ we de- 
rive in this way the leading high-frequency behaviour of 
the kernel, 



D-2 



D/2 



a u B Ry (8) 



where as is the Bohr radius, = 23/30 and ap = 8/15 
in D = 3, while ap = 11/16 and ap = 9/16 in D = 2. 
The longitudinal component of this result was obtained 



by Glick and LongE3 in D = 3 and by Holas and Singv 
in D = 2. In the low-frequency limit instead, perturba- 
tion theory gives an artificial divergence in Ref xc and a 
discontinuity in Im/ IC , and will not be pursued further. 

We obtain an approximate nonperturbative evaluation 
of (0) by decoupling its RHS within an RPA-like scheme, 
which in the frequency domain gives 



Im((AB;CD)) o 



ImU-C))^ 



(9) 



xlm ((B: £>»„_„, + Im ((A; £>»„, Im ((B: C)) w _ u ,] . 

The functions in the RHS of Eq. ^ are understood to 
be exact response functions of the interacting electron 
gas. This scheme clearly neglects exchange processes in 
the final state, which are known to reduce the spectral 
strength by a factor of 2 at high frequency in perturba- 
tive treatments (see Appendix O , but are ineffective at 
low frequency. This can be physically understood as fol- 
lows. A two-pair excitation at long wavelength involves 
the creation of holes with momenta p and p' inside the 
Fermi sphere and of electrons with momenta p + q and 
p' q outside the Fermi sphere. If the excitation energy 
lo is small compared to the Fermi energy Sp, then nec- 
essarily |q| <C kp and, since on average |p — p'| ~ kp, 
each electron is substantially closer (in k-space) to "its" 
hole than to the other one: exchange processes are thus 
suppressed in this case (see Fig. ^.a) by a factor v^ p /v^ 
which vanishes for q — ► (i.e. ui — ► 0).E3 Conversely if 
uj ^> Ep one has |q| ^> kp and therefore for parallel spins 
the strengths of direct and exchange processes are equal 
and opposite (see Fig. pj.b). 

On the basis of the above argument, which can be made 
quantitative at a perturbative level, we phenomenologi- 
cally incorporate exchange effects by inclusion of a factor 



9x(u) 



1 + 0.5^/2^ 
1 + uj/2e F 



(10) 



which interpolates between 1 at low u> and 0.5 at high 
uj on the energy scale (2ep) of exchange processes. Our 



final expression for the imaginary part of the kernel thus 



j fij,Ti \ , s r du' f d D q 



L\2 



2 



aL.T-^ImXL{q,u') + b L , T -t-ImxT{q,uj') 

uj z ' UJ Z 



x t —Im X L(q, w - uj'), 

(UJ — uj'Y 



(11) 



with cll.t as defined below Eq. (^), = 8/15 and 
b T = 2/5 in D = 3 and b L = b T = 1/2 in D = 2. 
The expression for the longitudinal component is equi#-. 
alent to that obtained in 3D by Hasegawa and WatabeEil 
by diagrammatic mea^ts and similar to that derived in 
2D by Neilson et a/J23 (see the discussion in Ref. p0| ). 
This result can be physically understood as representing 
the spectral density of two-pair excitations, which are 
present at any frequency in the long-wavelength region 
of the spectrum and provide a mechanism for the decay 
of the plasmon outside of the single-pair continuum. The 
resulting linewidth is Tk/ujk = -Im/^wt)/^. 

Considering the low-frequency behaviour of 
Im XL,T(k,uj) it can be shown that at low uj Eq. ( pT| ) 
is linear in uj, and that to this order the first term in the 
square brackets does not contribute. The coefficients of 
this linear behaviour are relatedE-3 to the bulk and shear 
viscosities of the electron gas, Q and rj respectively, via 



c = 



-n 2 lim 

ui— >0 



and 



V = 



-n lim 

UJ^O 



D 



Im /J c (w) 



(12) 



(13) 



The significance of— C and 77 is the same as in classi- 
cal hydrodynamics.^ Both viscosities vanish in the ideal 
Fermi gas, as well as in the RPA and in static-local-field 
approximations to the interacting electron gas. Since 
bh/br = 2(£) — 1)/D, also within the present model the 
bulk viscosity £ vanishes identically. Numerical results 
for the shear viscosity 77 will be given below. 

We now discuss the real part of the kernels /^(uj). 
Once the imaginary part has been evaluated, the real one 
can be obtained by means of the Kramers-Kronig relation 



Re/ x V H = ti? (00 



(14) 



where f^ T (oo) denotes the (real) high-frequency limit 
of the kernels. This quantity corresponds to the long- 
wavelength value of tha— leading spectral moment be- 
yond the /-sum ruleJ^J'O which is the first term on 
the RHS of Eq. (§). With the definition M L . T {k) = 
\\va u ^ OQ uj 2 XL,T{k,uj), we have 
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M L (k) 



nk 
2m 2 
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2m + D 
(k • q) 2 



fc 4 



(fee) 

OS(|q + k|) -£(?)) 



(15) 



and 
M T (fe) 



nfc 2 
2m 5 



2 (k • q) 2 



fe 4 



(5(|q + k|)-5( ? )) 



(16) 



where 5(g) is the static structure factor and (fee) denotes 
the true kinetic energy per particle. Expansion at long 
wavelengths gives the required high-w limit of /k C ' t (w), 



d L , T ((ke)-(ke)°) + e L , T ( P e)] (17) 



with d L = 4, c? T = 4/3, e L = 8/15 and e T = -4/15 in 
D = 3. and dj, = 6, c?t = 2, = 5/4 and er = —1/4 in 
D = 2. The average kinetic and potential energies (fee) 
and (pe) can be obtained from the Monte Carlo equation 
of state, and (fee) denotes the ideal-gas value. The 
resulting values of f^ T (oo), given in Table | (D = 3) and 
Table || (D = 2), allow to evaluate numerically Eq. Jl4| ) 
and to obtain the real part of the kernels at any fre- 
quency. !_. 

The low- frequency limit of f X c T {w) is relatedO to the 
elastic moduli K and /i via 



K xn = n 



/«c(0)-2^i/J c (0) 



and 



2 /J c (o). 



(18) 



(19) 



The significance of K and \x is the same as in classi- 
cal elasticity,!^ and as usual the suffix xc indicates that 
the ideal-gas contribution has been subtracted. Our 
results for K xc are in good agreement with the accu- 
rate values -obtained from the Monte Carlo xc energy 
per pgsticleB-y ef c c (n) via K™ c = n 2 d 2 ef c c / dn 2 (see 
belowE-a). We finally note that K xc is also related to the 
long-wavelength limit-of the static f xc (k,0) via the com- 
pressibility sum rule,E3 



A" 



lim lim f xc (k, 



k—>Q u— >o 



U3) 



(20) 



III. NUMERICAL RESULTS WITHIN THE RPA 

This Section presents results that we have obtained 
for f xc from numerical integration of Eqs. (11) and Jl4| ) 
using RPA response functions, which are given by 



vt T . (21) 



xJ^rH&i w ) + n/m %° T (fc, w) + n/m lu 2 k 

In the high-fre quen cy limit the first term in the square 
brackets of Eq.(|ll]) dominates over the second one, im- 
plying Im/J c = (cit / a L)Im f x c ; conversely, for low uj the 
first term is negligible and Im/J c = (6y/bi)Im/^ c . In- 
deed the longitudinal and transverse spectra are very 
similar, and the transverse one is rather accurately re- 
produced at all frequencies by setting Im/J c (w) ~ 0.72 ■ 
Im/^w) in 3D and Im/J c (w) ~ 0.85 • Im/£(w) in 2D, 
the proportionality factor being close both to ar/aL and 
to br/bL- For the real part of the kernels there is an 
additional shift due to the different values for u) = oo. 

Figure || reports the results for the imaginary part of 
fxc T in 3D at r s = 3 as functions of lo/lu p i [in 3D r s 
is defined as (47rna^/3) _1 / 3 and the plasma frequency 
LUpi as (47re 2 n/m) 1 / 2 ]. Bath, our results and the Gross- 
Kohn (GK) interpolatiorotHl for Im/'/ reproduce the 
high frequency behaviour given by Eq. (|^) and are linear 
at low frequency. As already remarked the bulk viscosity 
£ vanishes identically within the present model, but the 
shear viscosity rj is finite; numerical results for the latter 
are given in Table ffl. From Fig. || one can see that our 
estimate for 77 is significantly smaller than that of GK. 
In the intermediate frequency region our curves exhibit 
a sharp threshold at twice the plasma frequency, which 
is due to the onset of two-plasmon processes. This is the 
most remarkable difference from the GK interpolation. 

The corresponding real parts of the kernels are shown 
in Fig. ||. The two-plasmon threshold in the spectrum 
generates a pronounced minimum at lo — 2lu p i in the 
real part, which is absent in the smooth GK interpola- 
tion. The /4(0) value was obtained by GK assuming 
/|(0) = K xc /n 2 , i.e. /J c (0) = fee = [see Eqs. (|D| 
[n])]; our curves instead are consistent with Eq. ( |l8| ) with 
a nonvanishing f xc {0)- Table | reports our results for 
the elastic moduli K xc and fi xc , obtained from Re/ JC via 
Eqs. ( p^[jl^ ). They cannot be expected to be very pre- 
cise, being obtained through integration over the entire 
spectrum. With this caveat, we note that our numerical 
estimates for K xc agree with the accurate Monte Carlo 
results K^F , also given in the same Table, within 5%. 
We note that the xc contribution to the shear modu- 
lus [i xc becomes negative at low density; the total shear 
modulus [i = (i xc + %nep , however, remains positive. 

In order to facilitate use of these data in practical TD- 
DFT computations, an analytical fit has been given in 
Ref. |3|. The constraint f xc (0) = K xc /n 2 , which was 
imposed in the fits following GK, can be removed simply 
by setting (3=1. 

We now turn to the two-dimensional system. Figures ^ 
and H report our results for the imaginary and real parts 



of f xl r at r s = 3 [in 2D r s is defined as 



una. 



-l/2]_ 

The main difference with respect to the 3D results is the 
absence of the sharp two-plasmon threshold in the two- 
pair excitation spectrum (i.e. in Imf xc ), due to the fact 
that the plasmon dispersion relation vanishes in 2D at 
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long wavelength. Correspondingly the minimum in the 
real part of the kernel is much broader. The figures also 
compare our result for /,^_with the interpolation scheme 
of Holas and Singwi (HS) jEl which is the 2D extension of 
the GK interpolation. Both curves reproduce the asymp- 
totic limit (|l7|) as well as the high-frequency behaviour 
of Eq. (||), and both imaginary parts are linear in uj at 
low frequency. Also in the 2D case the minimum in the 
real part is absent in the GK/HS interpolation scheme. 

Table || reports the resulting values of the shear vis- 
cosity 77 and of the elastic moduli K xc and /j. xc , obtained 
as in the 3D case. In 2D the agreement between the bulk 
modulus K xc obtained from our data on f xl i T and the 
Monte Carlo value K^F is not as good as in 3D, but still 
better than 10% at all values of r s that we have consid- 
ered. In contrast to 3D, at low density the total shear 
modulus /1 = /i xc + ^uef, turns out to be negative. How- 
ever, the observed disagreement between K xc and 
prevents us from drawing a conclusion about the presence 
of an instability. 

To facilitate application of these results in actual TD- 
DFT computations, such as the one of Ref. [l^, we provide 
below a simple analytical interpolation. The imaginary 
part is reproduced by 



-gx(u) 



c\uj + c 2 uj 2 + C3UJ 3 + 2chs^ 5 

Co + C4LU 4 + UJ 6 



(22) 



where uj is in Ry, f xc in units of Ry/n and chs = Hir/Sr 2 
is the coefficient of the leading high-frequency behaviour 
from Eq. (^|). The remaining parameters, which we ob- 
tained by a least-squares fit to the numerical calculation, 
are reported in Table [II; the transverse component can 
be approximated by Im/J c = 0.85 Im/r c . The quality of 
a typical fit is shown in Figures || and Values at inter- 
mediate r s are best obtained by interpolation of f xc (uj) 
at the same uj (in Ry); for r s > 3 also the simpler scheme 
of interpolating the fitted coefficients cq, . ..,04 is viable. 



The real part can be obtained from Eq. (14), which in 
this case can be integrated analytically; the resulting ex- 
pression is quite long and we do not report it. 



IV. EXCHANGE-CORRELATION EFFECTS 

ON f xc 

This Section is aimed at assessing the validity of the 
results presented above. In the first part we study the ef- 
fect of correlation in the treatment of each pair, adopting 
more refined response functions in the RHS of Eq. (11); 
in the final part we discuss corrections which go beyond 
the present two-pair model. 

We have introduced the effect of correlation in the 
treatment of each pair in Eq. ([ll]) by means of two of 
the most successful static-local-field approximations, the 
SingwipTosi-Land-SjolandeiO (STLS) and the Vashishta- 
SingwiE^I (VS). In the 3D case, both schemes predict neg- 
ative plasmon dispersion at large r s (r s > 5 in STLS and 



r s > 9 in VS), in qualitative agreement with experiment .0 
The VS scheme embodies the compressibility sum rule on 
the static response, and is therefore more reliable in the 
study of static phenomena. Since these schemes only in- 
volve longitudinal currents, the transverse response func- 
tion is still evaluated at an RPA level. 

Figures |and| compare the RPA results in 3D with 
those obtained with STLS and VS. At r s = 1 correlation 
only gives minor corrections, but at larger r s it leads 
to a divergence, caused by the appearance of negative 
plasmon dispersion. Figure |^ compares the results ob- 
tained with STLS at various densities. The minimum 
at intermediate frequency gets more pronounced with in- 
creasing r s , and for r s > 5 becomes a divergence of the 
form 6(uj - 2uj min )(uj - 2uj mm )~ 1/2 , where uj mm is the 
minimum energy of the collective mode. 

Introduction of the local field correction significantly 
increases the shear viscosity 77, and leads to large negative 
values of the shear modulus fi xc (e.g. with VS we get 
fj, xc = —0.03 at r s — 5 and fi xc = —0.1 at r s = 10, in 
units of 2uj p in). The relation K xc — K^F is satisfied 
with the same accuracy as in the RPA case. 

Figure ^| displays the dynamical structure factor 
S(k,uj) — — (2k 2 /nuj 2 )ImxL(k, uj) which is relevant to in- 
elastic scattering experiments. The threshold behaviour 
at frequency 2uj p i is a clear-cut signature of the present 
results for f xc . The more pronounced singularity present 
in the STLS curves originates from the negative plasmon 
dispersion. 

We also investigated the role of exchange and correla- 
tion in 2D r _Lising the STLS model as generalized in 2D 
by JonsonO The qualitative behaviour turns out to be 
similar. Figure [l0| reports our results for the imaginary 
part of f xc as obtained from RPA and STLS calculations. 
As in 3D, in STLS the plasmon energy at intermediate 
wavevector is lower than in RPA, and correspondingly 
we note a significant increase in the depth of the mini- 
mum at intermediate frequency in Imf xc . This increase, 
unlike the 3D case, is significant also at rather high den- 
sity (see the inset) as a consequence of the enhancement 
of correlation effects in 2D systems. On the other hand, 
in 2D no divergence appears in the two-pair spectrum, 
since no minimum at finite k is present in the plasmon 
dispersion. The same is true for the corresponding real 
parts, which are shown in Fig. pT[ 

Having discussed at length the role of two-pair pro- 
cesses in the long-wavelength spectrum of the uniform 
electron gas, we now turn to a qualitative discussion of 
other effects which have been neglected in the present 
treatment. Whereas a detailed quantitative calculation is 
cumbersome, the qualitative role of multi-pair processes 
beyond two-pairs can be easily grasped by straightfor- 
ward extension of the present treatment, in the spirit of 
an expansion in the number of pairs involved. A per- 
turbative analysis shows that in the high-frequency limit 
these processes are of higher order in 1/uj, and therefore 
do not contribute to the asymptotic behaviours of Eq. 
(0). In the intermediate-frequency regime one can ex- 
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pect the appearance of a n-plasmon threshold effect in 
the n-pair channel. Whereas at present we do not have a 
reliable quantitative estimate of the contribution of such 
processes to Imf xc , we believe that their overall spectral 
strength will be a minor correction to the present results. 
This can be inferred from the good agreement between 
the two-pair result for K xc and that from the Monte 
Carlo data, the difference being a quantitative measure 
of the integrated spectral strength of higher-order pro- 
cesses. 

A different class of effects which could modify the 
present results are improved response functions in the 
two- pair channel, i.e. in Eq. (pi]). The main qualita- 
tive feature which is absent from all response functions 
that we have considered is plasmon damping. This will 
broaden the sharp feature present in 3D at the two- 
plasmon threshold. When the plasmon dispersion is posi- 
tive (local-field-corrected results at small r s and RPA) we 
expect small corrections, since the threshold behaviour is 
driven by the long- wavelength plasmon, whose linewidth 
vanishes as k 2 for small k. In the case of negative plasmon 
dispersion, the collective excitation with minimum en- 
ergy has nonzero wavevector and therefore non- vanishing 
linewidth due to decay into multiple particle-hole pairs. 
This indicats that the divergence in the two-pair spec- 
trum found in the local-field-corrected results at large r s 
(see Figures ^ and ||) is probably an artifact and would 
be replaced by a smooth peak if plasmon damping were 
included. 



V. SUMMARY 

In this work we have given an extensive discussion of 
the exchange-correlation kernels f^ T (u!), both in 2 and 
in 3 spatial dimensions. We have presented an exact ex- 
pression for the kernels in terms of four-point response 
functions, and evaluated it numerically within a non- 
perturbative approximate decoupling scheme, which ac- 
counts for two-pair processes. Our numerical results are 
qualitatively different from previously known interpola- 
tions. In 3D we predict a threshold behaviour, which can 
be understood as due to a simple phase-space effect, i.e. 
the opening of the two-plasmon channel in the two-pair 
spectrum. In 2D the same mechanism leads to a broad 
feature in the spectrum. We have also studied the in- 
fluence of static-local-field corrections on our results and 
found even more marked deviations from previous theo- 
ries. 

We have obtained good agreement with all known 
asymptotic behaviours, including the newly-derived high- 
frequency limit of the transverse part and the Monte 
Carlo results for the bulk modulus. Estimates for 
the shear modulus and viscosity have also been given. 
Our results, which have been fitted to simple analyti- 
cal formulas, provide the entire input necessary for non- 
adiabatic TD-DFT computations. 
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APPENDIX A: TRANSVERSE RESPONSE 
FUNCTION FOR THE NON-INTERACTING 
FERMI GAS IN 2D AND 3D 

The current-current response function for the non- 
interacting Fermi gas is given by 



$(k, W )=£ 



'q+k/2 



Vk/2 



m 2 lu ~ (e q - k /2 ~ Eq+k/2) 



:, (Ai) 



where = 28(kp — |q|) is the Fermi distribution func- 
tion and e q — q 2 /2m the free particle energy. It is suffi- 
cient to evaluate the transverse component of 



Aij (k, to) 



d D q F (q+ik) ,: (q+Ik) 



{2irf q "< 
x6(uj + e q - e q+k ) 



(A2) 



since Imx^k, u) = (k, u>) — (k, — ui). Angular- 
integration leads to 



A T (k,u) = -■ 



1 



-e fc |/fe 



dq 1D q D 6{k F - q) (A3) 



where 72 = \/l — x 2 , 73 = (1 — x 2 )/4 and x = m\uj — 
Ek\/qk. It is convenient to work with the reduced vari- 
ables z = k/2kp and u — Lum/kkp. After integration we 
have 



Tm]&(k,w) = l D — (B+-B-) 



(A4) 



la = 



with B± = 0(1- I u± z |) l-(u±z) 

3tt/32 and l 2 = 1/3. 

The imaginary part of the ideal-gas transverse response 
function vanishes when \ \tu\ — 2kkF/m\ > Ek, as the lon- 
gitudinal one. As a check of the present derivation we 
evaluate the third frequency moment sum rule, which 
leads to 

2 r°° n 

M$(k) = / wlmxr(ft,u)(iu = h D — e k e F (A5) 

71" Jo m 



with hs = 4/5 and I12 = 1, in agreement with (|16|) eval- 
uated to zero order in the interaction potential vi , i.e. 
with = 0, (fee) = fep in 3D and (he) = \e F in 2D. 



G 



The real part is obtained from the imaginary one via 
the Kramers-Kronig relation, which gives 



3n 
Smz 
1 



z d + 3u 2 z z 

3 



with E. 



(3) 



[(z±- 



4 E 



l] 2 ln 



(3) 



1 

4 E 



(3) 



u±z — 1 



u±z+l 



in 3D, and 



Rex T (k,uj) 



Smz 



2z 3 + 6u 2 z -3z- E. 



(2) 



E 



(A6) 



(2) 



(A7) 



with E { 1 ] = [{z ± u) 2 - 1]30(| z±u\ -l)sgn(z ± u) in 
2D. 

We examine now some limiting behaviours. At long 
wavelength we get 



n „ x Af£(fc) 



9d- 



-F 2 F 2 
i tit j 



"T c fc c F ' 



(A8) 



where g 3 = 48/35, g<z = 1/2 and M^(fc) as defined in 
(A5). At zero frequency we have instead 



Xr(*»w = 0) 



n 


'5 


3 k 2 


m 


8 


32 /c^r 


3fc 


F 


^ fc 2 


8 fc 





2 

1 I In 



k - 2fc F 



+ 2k F 



(A9) 



in 3D, and 

Xt(*,w = 0) 



n 



64 



~3fc \4fc 2 



-1 6(k-2k F ) 



(A10) 



in 2D. In both cases the w — > 0, fc — * limit depends on 
the order in which the limits are taken, as in the longi- 
tudinal case. 



APPENDIX B: EVALUATION OF THE 
LOW-fc EXPANSION 

In this Appendix we evaluate to leading order in k 
the second term in the RHS of Eq. (0) , which obeys the 
equation of motion 

cu 2 Im «[j k) H]; [j_ k , H])) u = Im (([[j k , K + P], K + P]; 

[[j_ k ,K + P],K + P]» w , (Bl) 

where K and P denote the kinetic and potential terms 
of the Hamiltonian respectively. There are 16 terms in 
total, which correspond to the different combinations of 
K and P. Fortunately only 4 of them are relevant to our 



calculations, as is shown in the following. Every commu- 
tator with K gives a factor of fc, therefore terms contain- 
ing more than two K's do not contribute to leading order 
in fc. The commutator [[j k ,P],P] vanishes, since 



[jk,P] = [jk.J^PqP-q] 
q#0 



V ^ m 

q#0 



Pq+kP-q (B2) 



and two density operators p q commute. 
The remaining terms can be written as 



1 



Im(([[j k ,K],P] + [[j k ,P],K] ; 
[[j_ k ,K],P] + [[j_ k ,P],K])) u 



(B3) 



At this point we remark that the term P 

4 



V lv kPkP-k OI the Hamiltonian gives rise to a singu- 



(B4) 



lar contribution to the equation of motion for j k , 

fc 2 

[k • j k , P rcs ] = ufypk - — v% p 2k/ o-k. 

In the first term in the RHS the factor V^ 1 present in 
the others has been compensated by the factor po = N. 
Therefore this term has to be treated separately when 
transforming the summations into integrals. We found 
it more convenient to replace the Hamiltonian H with 
H = H — P ICS in the equation of motion for the re- 
sponse function, which allows to obtain the proper re- 
sponse function \- This can be easily understood from 
the diagrammatic definition of \ as the sum of all dia- 
grams that cannot be split by cutting a single interac- 
tion line: this necessarily carries a factor Wq =k , which 

has been excluded from H . Using H one also excludes 
all non-singular contributions having a factor fq_ k , as 



for example the second term in the RHS of (B4). These 
terms however contain a factor V~ x and are irrelevant in 
the thermodynamic limit. 

All considerations made previously hold with the new 
Hamiltonian, provided that we are computing the proper 
current-current response in Eq. (|^). The commutators in 
(B3) are given by 



q V / 

[[Jk, P] . K] = ^ £ [if q+t| (q + k)< - <q 



x(q + k) • j q+k p_ q . 



kP-q • 



(B5) 



where terms containing a single p k operator have been 
neglected for reasons explained below. From Eqs. (|q), 
(pi and p| we get 



Imxii(k,w) = -— L_^Im((j q+kP _ q ;j^ 



m 2y2^4 \\Jq-t-Kr--MJ Jq'-kP-q 

qq' 

xr li (q,k)r^'(q',-k)+ (fc 2 ) 



(B6) 
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where T is given by 



With a, 6, <7, h < kp and c, d,e, f > hp, we have 



r«(q,k) 



W |q+k| 
L („l 



qV-k'q'-^qk 



(B7) 



We remark that ( |B6| ) is valid for every Fourier- 
trasformable interparticle potential, and is not restricted 
to T = 0. 

The tensor T is of fi rst order in k and the product of 
the two tensors in (B6) is of second order, so that in the 
response function we can safely set k — 0, leading to 
Eq. (^) . This also explains why terms co ntai ning one pk 
have been neglected in the commutators (B5): po is con- 



stant in time and its commutator with j q /o_ q vanishes. 



APPENDIX C: EVALUATION OF THE 
HIGH-FREQUENCY LIMIT VIA 
PERTURB ATIVE EXPANSION 

Starting from Eq. (^) we evaluate the high frequency 
limit of fxc{u) to second order in perturbation theory. 
In this way we extend to the-transverse case the results 
obtained by Glick and Longc3 in 3D and by Holas and 
SingwiEZI in 2D on the asymptotic behaviour of f xc . The 
discussion will also clarify the nature of the exchange pro- 
cesses which we approximately included in the decoupling 
(H|) via the factor g x of Eq. ©. 

The four-point response function can be written as 



Im ((jqP-q;jq'P-q' 



(Pl + f) i (p 2 -4) i 



P1P2P3P4 



X ^5(w„ - W)(0|c+ C+ 3 C pi+q Cp3_ q | 



x(n\c + ,c + 

\ l^p 2 -q' p 



4+q' C P2 



= P4 |o) 



(CI) 



where |0) and \n) denote the exact eigenstates of the 
system, and the spin index is implicit. Since each T in 
(B6) contains a factor vjf, we can evaluate the four-point 
response function at zero order in perturbation theory, 
i.e. for the non-interacting electron gas, so that lu u q — 
(q.' 2 /m) + (p 4 - p 3 ) • q'/m. 

The product of expectation values in the previous 
equation is different from zero only when P\,P2tP3iPa < 
kp, which in the high-frequency limit implies p t <C q — 
q' ~ y/rrnj. These considerations allow us to simplify 
©as 



(0\ctc^c c c d cfc^c g c h \0) = {S ath S big - S a<g S bt h) 

x {5d,e5 c j - S c , e 6dj)- (C3) 

There are 4 terms overall, two of which are negative be- 
cause of the anticommutation rules and mix together field 
operators belonging to different density (or current den- 
sity) operators. These terms, which we call exchange 
terms, are neglected in the decoupling (|9|) and carry an 
overall factor — \ due to spin (see below) . The first term 
5a,h$b,g&d,e$c,f imposes the following constraints on the 
wave vectors and spin indices, 



Pi 

P2 



P4 CTi 
P3 0"2 



04 Pi 
CT3 P3 



q = P4 
q = P2 



q 0i 

q' cr 2 



04 , 
0"3 ■ 



(C4) 



The sum is done on pi,0i P2,0'2 J and brings a factor 
N 2 Sq,qi, with N the number of particles. The second 
term 5 a ,g6 bth 5d,e$cj gives 



Pi = P3 Cl 
P2 = P4 0"! 



03 
CT 4 



Pi 
P3 



q 
q 



P4 
P2 



(T 2 

a 2 



04 , 
0"3 ■ 



(C5) 



Unlike the first one, this second term fixes all the four 
it's, contributing with a factor — ^-£ q , q ' which is half 
the opposite of the first one. 

The remaining two terms are obtained in an analo- 
gous way, and their overall contribution is also ^-5, 
Putting all terms together we get 



q>-q 



Im 



l . •/' 
JqP-qi Jq'P— q' 



^ (dq,q' 

,2 



2 m 2 



q'q' ( ' x ( q 

X LO 

4 \m 



(C6) 



Substitution of the previous expression in Eq. (^|) and 
small k expansion lead to Eq. (||). 

An estimate of higher-order corrections can be ob- 
tained from additional applications of the equation of mo- 
tion for j q . Third-order perturbation theory treats one of 
the two pairs to first order and is irrelevant for large q and 
u>. To fourth order in Vu, one obtains instead three-pair 
contributions which contain two additional commutators 
with P in the last term of Eq. (j^) . Considering that the 
transferred momentum scales as y/U at large uj, these 
scale as and therefore do not affect the leading 

behaviour given by second-order perturbative expansion. 



Im 



JqP-qJq'P-q' 



q'q"' fq' 2 

7T- — d Ul 

4m z \ m 



E ( l C p"l C P3 C Pl+q C P3-q C p 2 -q' C p4+q' C P2 C P4l ) ( C2 ) 



P1P2P3P4 
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"/" P i 

FIG. 2. Imaginary part of f^ T {u>) in 3D at r a = 3 in units 
of 2ui v i/n, asriurictions of u/u) p i. The dashed line gives the 
Gross-Kohnu't-J interpolation scheme. 




1 2 3 4 5 

"/" P i 



FIG. 3. Real part of fb T (u) in 3D at r s = 3. Notations 
and units are as in Fig. ^. The dot on the left axis marks 
K™ c c /n 2 . The scale for the transverse component is on the 
right-hand axis. 



2 4 6 8 

w(Ry) 

FIG. 4. Imaginary part of /ic T (o;) in 2D at r B — 3 in units 
of Ry/n, as fptictions of u (in Ry). The dashed line gives the 
Holas-SingwiEj interpolation and the dotted lines are the fit 
discussed in the text. 
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FIG. 5. Real part of fxc T ( l ^ J ) m 2D at r B = 3 in units of 
Ry/n, as functions of ui (in Ry), the scale for the transverse 
component being on tha right-hand axis. The dashed line 
gives the Holas-SingwiEZI interpolation and the dotted lines 
show the fit discussed in the text. The dot on the left axis 
marks Kf^/n 2 . 
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FIG. 6. Imaginary part of f^ c (uj) in 3D in units of 2u p i/n, 
as a function of lo/lo p i at r s = 5 (main figure) and r s = 1 (in- 
set). We plot results obtained using RPA response functions 
(dashed curve), STLS response functions (full curve) and VS 
response functions (dotted curve). 




FIG. 7. Real part of /^(w) in 3D in units of 2oj p i/n, as 
a function of lo/lo p i at r s = 5 (main figure) and r s = 1 (in- 
set). We plot results obtained using RPA response functions 
(dashed curve), STLS response functions (full curve) and VS 
response functions (dotted curve). The dots mark K^P/n 2 . 
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FIG. 8. Imaginary part of f^ c {uj) in 3D, in units of 
2r^ 3 ^ 2 u> p i/n, computed with STLS response functions, at var- 
ious values of r s . 
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FIG. 9. Dynamic structure factor S(k, u>) at r s — 5 as a 
function of uj/uj p i shown in the flags on a semilogarithmic 
scale at various values of kr a an, as obtained from STLS (full 
curves) and RPA (dashed curves) calculations. 
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FIG. 10. Imaginary part of f^ c {uj) in 2D in units of Ry/n, 
as a function of u> (in Ry) at r s — 6 (main figure) and 
r s = 1 (inset). We plot results obtained using RPA response 
functions (dashed curves) and STLS response functions (full 
curves) . 
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FIG. 11. Real part of /^(w) in 2D in units of Ry/n, as a 
function of w (in Ry) at r s = 6 (main figure) and r s = 1 (in- 
set). We plot results obtained using RPA response functions 
(dashed curve) and STLS response functions (full curves). 
The dots mark Kf c c /n 2 . 
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TABLE I. Exact high-w limits of f^ T (u) in 3D and bulk modulus K^P obtained from the 
Monte Carlo equation of state; shear viscosity rj and shear and bulk moduli (ji xc and K xc ) obtained 
from the RPA treatment of two pair processes. f xc is in units of 2ui p i/n, rj in units of n, K xc 
and /j, xc in units of 2u) v in. Values in atomic units (as — rn = e 2 = 1) can be obtained from 
r] AU = 3/(47rrJ)r; tab , (K,fi) AU = 3 3/2 r7 9/2 /(2n)(K, ji) tah , where tab denotes the tabulated values. 
The last reported decimal figure is likely to be affected by numerical inaccuracies. 



r s 


t-MC 
*^ x c 




/Jc(oo) 


n 




K 


0.5 


-0.04246 


-0.01794 


0.0177 


0.0029 


0.0065 


-0.0425 


1 


-0.0611 


-0.0216 


0.0284 


0.0062 


0.0064 


-0.0612 


2 


-0.0891 


-0.0252 


0.0457 


0.012 


0.0052 


-0.0896 


3 


-0.1119 


-0.0280 


0.0600 


0.017 


0.0037 


-0.1128 


4 


-0.1320 


-0.0308 


0.0724 


0.021 


0.0020 


-0.1335 


5 


-0.1503 


-0.0338 


0.0835 


0.024 


0.0002 


-0.1525 


6 


-0.1674 


-0.0370 


0.0935 


0.027 


-0.0018 


-0.1702 


10 


-0.2276 


-0.0518 


0.1267 


0.034 


-0.010 


-0.233 


15 


-0.2917 


-0.0725 


0.1587 


0.040 


-0.023 


-0.301 


20 


-0.3483 


-0.0939 


0.1847 


0.044 


-0.036 


-0.361 


TABLE II. Exact high- 


-lo limits of f xc 


T {k,uj) in 2D and bulk modulus K^F obtained from the 


Monte Carlo equation of state; shear viscosity rj and shear and bulk mo 


duli (ji xc and K xc ) obtained 


from the RPA treatment of two-pair processes. f xc is in units of Ry/n, rj in units 


of n, K xc 


and ji xc 


in units of Ry ■ 


n. Values in 


atomic units can 


be obtained from rj AU = l/^r 2 )^ 15 , 


(K,ji) AU = l/(27rr 2 )(A» tab , where tab 


denotes the tabulated values 


. The last reported decimal 


figure is 


likely to be affected by numerical inaccuracies. 








r s 


k mc 


f X c(oo) 


/lc(°0) 


V 


jl X c 


Kxc 


1 


-0.9360 


-0.5499 


0.3372 


0.018 


-0.064 


-0.959 


2 


-0.4912 


-0.2750 


0.1916 


0.029 


-0.064 


-0.514 


3 


-0.3413 


-0.1933 


0.1330 


0.035 


-0.058 


-0.363 


4 


-0.2649 


-0.1535 


0.1010 


0.040 


-0.054 


-0.285 


5 


-0.2180 


-0.1294 


0.0810 


0.043 


-0.050 


-0.236 


6 


-0.1860 


-0.1128 


0.067 


0.045 


-0.047 


-0.203 


10 


-0.1191 


-0.0768 


0.0395 


0.050 


-0.039 


-0.132 


15 


-0.0833 


-0.0560 


0.0257 


0.052 


-0.033 


-0.094 


20 


-0.0645 


-0.0445 


0.0189 


0.054 


-0.028 


-0.073 



TABLE III. Interpolation parameters according to 
Eq. ([22J), in 2D. The last reported decimal figure is likely 
to be affected by numerical inaccuracies. 



r 3 


1n -3 5/2 
10 C JV 


cir 2 


C2 


C3 


c 4 


1 


62.7 


1.10 


9.94 


37.4 


6.84 


2 


2.90 


59.2 


-1.74 


4.62 


7.70 


3 


1.73 


34.9 


-3.79 


12.5 


15.1 


4 


1.44 


28.6 


-3.45 


15.8 


30.3 


5 


1.18 


22.5 


-2.64 


16.3 


47.9 


6 


0.943 


17.3 


-2.04 


15.8 


66.2 


10 


0.458 


7.23 


-1.0 


13.1 


150 


15 


0.250 


3.42 


-0.524 


10.9 


276 


20 


0.160 


1.96 


-0.313 


9.52 


430 



13 



